wd<- "/Users/chengg/Google Drive/EPICON/Mycobiome/Fungal ITS/statistic/Total.fungi"
setwd(wd)
library(vegan)
library(ecodist)
rm(list = ls())
rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ghost.2019.03.17.Rdata")
fung<-fung.rar #rarefied fungal communities
env<-env # meta data
plot.gd<-function (habitat, j, k, z) #create a 'plot.gd" function
#habitat: sample types (Root, Rhizosphere, Soil, Leaf)
# k: the y axis value to print mantel result, 0.1 to 0.9
# the abbreviation of the habitat (Leaf, Root, Rhizo, Soil)
{
fung1<-fung[env$Habitat==habitat,]
env1<-env[env$Habitat==habitat,]
par(mfrow=c(3,6),mar=c(2, 2, 0.5, 0.5))
for (i in c(j:17))
{
fungx<-fung1[ env1$TP== i,]
envx<-env1[ env1$TP== i,]
fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
gd<-distance(data.frame(envx$X, envx$Y))
td<-distance(envx$TP)
color=rgb(0,0,0,alpha=0.5)
plot(jitter(gd),fung.dist,xlab=" ",ylab=" ",,xlim=c(0,80),ylim=c(0,1), col=color, cex=2)
text(40, 0.9, sprintf("%s-week%s", z, i), col="red", cex =2.5, font =2)
ma<-mantel(fung.dist~gd, nperm = 999)
text(40, k, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
}
}
plot.gd("Leaf", 1, 0.7, "Leaf")

plot.gd("Root", 1, 0.1, "Root")

plot.gd("Rhizosphere", 1, 0.1, "Rhizo")

plot.gd("Soil", 0, 0.1, "Soil")

rm(list = ls())
rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ghost.2019.03.17.Rdata")
fung<-fung.rar #rarefied fungal communities
env<-env # meta data
plot.td<-function (habitat, j, k, z) #create a 'plot.td" function
#habitat: sample types (Root, Rhizosphere, Soil, Leaf)
# k: the y axis value to print mantel result, 0.1 to 0.9
# the abbreviation of the habitat (Leaf, Root, Rhizo, Soil)
{
fung1<-fung[env$Habitat==habitat,]
env1<-env[env$Habitat==habitat,]
par(mfrow=c(3,6),mar=c(2, 2, 0.5, 0.5))
for (plot in c(1:6,11:22))
{
fungx<-fung1[ env1$plot== plot,]
envx<-env1[ env1$plot== plot,]
fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
gd<-distance(data.frame(envx$X, envx$Y))
td<-distance(envx$TP)
color=rgb(0,0,0,alpha=0.5)
plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1), col=color, cex=2)
text(9, 0.9, sprintf("%s-plot%s", z, plot), col="red", cex =2.5, font =2)
ma<-mantel(fung.dist~td, nperm = 999)
text(9, k, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
}
}
plot.td("Leaf", 1, 0.1, "Leaf")

plot.td("Root", 1, 0.1, "Root")

plot.td("Rhizosphere", 1, 0.1, "Rhizo")

plot.td("Soil", 0, 0.1, "Soil")

Temporal distance
rm(list = ls())
rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ghost.2019.03.17.Rdata")
fung<-fung.rar[env$TP>0,] #rarefied fungal communities
env<-env[env$TP>0,] # meta data
env$Treatment2<-factor(env$Treatment1, labels = c("CON", "POST", "PRE")) #"CON", "POST", "PRE" are the abbre of Control, Pre_flowering_drought, Post_flowering_drought
par(mfrow=c(3,4),mar=c(2, 2, 0.5, 0.5))
for (treat in c("CON", "PRE", "POST")) {
fung1<-fung[env$Treatment2== treat,]
env1<-env[env$Treatment2== treat,]
for (habitat in c("Leaf","Root", "Rhizosphere", "Soil")){
fungx<-fung1[env1$Habitat==habitat ,]
envx<-env1[env1$Habitat==habitat ,]
fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
gd<-distance(data.frame(envx$X, envx$Y))
td<-distance(envx$TP)
color=rgb(0,0,0,alpha=0.5)
plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1), col=color, cex=0.5)
text(9, 0.95, sprintf("%s-%s", habitat, treat), col="red", cex =2.5, font =2)
ma<-mantel(fung.dist~td, nperm = 999)
text(9, 0.05, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
}
}

par(mfrow=c(1,4),mar=c(2, 2, 0.5, 0.5))
for (treat in c("CON")) {
fung1<-fung[env$Treatment2== treat,]
env1<-env[env$Treatment2== treat,]
for (habitat in c("Leaf","Root", "Rhizosphere", "Soil")){
fungx<-fung1[env1$Habitat==habitat ,]
envx<-env1[env1$Habitat==habitat ,]
fung.dist<-vegdist(decostand(fungx,"hellinger"),"bray")
gd<-distance(data.frame(envx$X, envx$Y))
td<-distance(envx$TP)
color=rgb(0,0,0,alpha=0.5)
plot(jitter(td),fung.dist,xlab=" ",ylab=" ",xlim=c(0,18),ylim=c(0,1), col=color, cex=0.1)
text(9, 0.95, sprintf("%s", habitat), col="red", cex =2.5, font =2)
ma<-mantel(fung.dist~td, nperm = 999)
text(9, 0.05, sprintf("R=%.3f, P=%.3f", ma[1], ma[2]), col="blue", cex =2, font =2)
}
}
